---
title: "Brumberg et al: Global Analysis of Constraints to Natural Climate Solution Implementation"
author: "Hilary Brumberg, Margaret Hegwood, Waverly Eichhorst, Anna LoPresti, James T. Erbaugh, Timm Kroeger"
date: "2025-05-22"
output: html_document
---

#Data Preparation

##Load packages
```{r}
library(readxl)
library(dplyr)
library(tidyr)
library(rnaturalearth)
library(countrycode)
library(tidyverse)
library(ggplot2)
library(ggthemes)
library(scales)
library(rcartocolor)
library(vegan)
library(pheatmap)
library(igraph)
library(ggraph)
library(tidygraph)
library(forcats)
library(grid)
library(gridExtra)
library(patchwork)
library(reshape2)
library(ggalluvial)

```

##Load and clean data
```{r}

#Read in data

rawDF <- read.csv(file.path(getwd(), "DatasetSI1.csv"))

placeCodes <- read.csv(file.path(getwd(), "countrycodes.csv"))

constraintCodes <- read.csv(file.path(getwd(), "ConstraintDefinitions.csv"))

Naturebase <- read.csv(file.path(getwd(), "Naturebase_clean.csv"))


#Clean data

placeCodes <- placeCodes[,-7]

df <- rawDF %>% filter(secondscreen == 1, reviewscreen == 0, yearpublication >= 2020) %>% dplyr::select(constraintID, countryconstraint, constraint, constraintpathway, solutions)

df <- df %>% mutate(countryconstraint = strsplit(as.character(countryconstraint), ", ")) %>%  unnest(countryconstraint)

df <- df %>% mutate(constraintpathway = strsplit(as.character(constraintpathway), ", ")) %>%
  unnest(constraintpathway)

colnames(df) <- c("constraintID", "ISOCode", "constraint", "pathway", "solutions")

df <- merge(df, placeCodes , by = "ISOCode")

constraintCodes <- head(constraintCodes, 39)

constraintCodes <- constraintCodes %>%
  mutate(constraintabb = str_trim(constraintabb))

df <- merge(df,constraintCodes, by = "constraint", all = TRUE)


df <- df %>% mutate(pathway = ifelse(pathway %in% c("GrR", "AGC"), "AGC/GrR", pathway)) #Group grassland NCS
df <- df %>% mutate(pathway = ifelse(pathway %in% c("PeR", "APC"), "APC/PeR", pathway)) #Group peatland NCS
df <- df %>% mutate(pathway = ifelse(pathway %in% c("CWR", "AWC"), "AWC/CWR", pathway)) #Group coastal NCS
df <- df %>% filter(pathway != "FPM", pathway != "RWF")

df <- df %>%
  mutate(paperID = sub("_.*", "", constraintID))

```



##Set colors
```{r}


#NCS Pathways
pathwayColors <- c("#AA4499", #Avoided Forest Conversion (AFC)
                      "#E6C20E", #Avoided Grassland Conversion (AGC)
                      "#44AA99", #Agroforestry (AgFo)
                      "#AF47E8", #Avoided Peatland Conversion (APC)
                      "#392599", #Avoided Wetland Conversion (AWC)
                      "#09712B", #Climate-Smart Forestry (CSF)
                      "#E8790D") #Reforestation (RFo)

#SDG Regions
sdgColors <- c("#D81B60", #Central and Southern Asia
               "#1E88E5", #Eastern and South-Eastern Asia
               "#D87A00", #Europe and Northern America
               "#7B037D", #Latin America and the Caribbean
               "#02940E", #North Africa and Western Asia
               "#1EA99B", #Oceania
               "#DAC200") #Sub-Saharan Africa

#Income Groups
incomeColors <- c("#9500FF", #high income
                  "#ffd166", #low income
                  "#06d6a0", #lower-middle
                  "#118ab2") #upper-middle

#Constraint Colors
constraint_colors <- c(
  #Economic
  "Lack of credit" = "seagreen1",                              
  "Lack of funding" = "green4",                          
  "Lack of insurance" = "darkgreen", 
  "Inadequate markets for ecosystem services" = "olivedrab2", 
  "Inadequate markets for outputs"  = "green",    
  "Inadequate prices for ecosystem services" = "olivedrab4",   
  "Inadequate prices for outputs" = "#4daf4a",  
  "Time lag in benefits"  =  "cadetblue3" ,                    
  "Tradeoffs with agriculture" = "mediumseagreen",    
  
  #Governments & Organizations
  "Corruption or lack of transparency"  = "yellow4",       
  "Lack of administrative capacity" = "lightgoldenrod",             
  "Lack of enforcement"  = "yellow",                      
  "Lack of coordination among organizations" = "gold2" ,
  "Lack of political will" = "lightgoldenrod4" ,                   
  "Violent conflict or civil unrest" = "khaki1" , 
  
  #Knowledge
  "Insufficient info about yields or profits" = "mediumpurple", 
  "Lack of info about co-benefits"   = "magenta4",           
  "Lack of info on how to design or manage" = "purple4",   
  "Inadequate planning or management" = "mediumorchid4" ,         
  "Lack of or insufficient technical support" = "mediumorchid" ,
  "Limited land manager capability"  = "purple1" ,          
  
  #Material Inputs           
  "Lack of input materials and infrastructure" = "mediumblue",
  "Lack of labor" = "skyblue",                              
  "Lack of water/water distribution networks" = "#377eb8", 
  "Unsuitable land" = "blue4",  
  
  #Policies & Rules
  "Ineffective laws, policies, or regulations" = "violetred",
  "Insecure land tenure" = "violetred4",                      
  "Unsecure or uncertain benefits" = "indianred",           
  "Incentives for non-NCS"   =  "pink",  
  "Lack of laws, policies, or regulations" =  "#e41a1c",  
 
  #Social & Behavioral
   "Concerns over negative equity impacts" =  "peru",     
  "Human-wildlife conflict" = "darkgoldenrod4",                   
  "Interpersonal conflict" = "lightsalmon4",    
  "Lack of opportunity to participate" = "#a65628",        
  "Limited social learning/exchange networks" =  "wheat4", 
  "Preferences for non-NCS land uses" = "tan",          
  "Disinterest or skepticism" = "brown",  
  
  #Negative Side Effects
  "Tradeoffs with other land uses" = "orange",            
  "Other negative side effects" = "orangered",  
  
  "Tie" = "grey40",
  "NA" = "white"
)

```


#Main Text Figures 

##Figure 1: 8 panel map of number of panels 
```{r}
map <- ne_countries(scale = "medium" ,returnclass = "sf")

fig1 <- df[, c("ISOCode", "constraintID", "pathway")]
fig1 <- fig1[order(fig1$constraintID),]
fig1 <- unique(fig1)

fig1a <- unique(fig1[,c("ISOCode", "constraintID")])
fig1a <- as.data.frame(table(fig1a$ISOCode))
colnames(fig1a) <- c("iso_a3", "freq")

fig1a <- merge(map, fig1a, by = "iso_a3", all = TRUE)

#ALL PAPERS
p1 <- fig1a %>% ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "lightgrey", 
                      high = "black", 
                      na.value = "white", 
                      limits = c(1,max(fig1a$freq, na.rm = TRUE))) + 
  ggtitle("A. All Pathways")

#Constraints by pathway by country
fig1b <- as.data.frame(table(fig1$pathway, fig1$ISOCode))
colnames(fig1b) <- c("pathway", "iso_a3", "freq")
fig1b$iso_a3 <- as.character(fig1b$iso_a3)

#AGROFORESTRY
agroforestry <- fig1b %>% filter(pathway == "AgFo")
agroforestry <- merge(map, agroforestry, by = "iso_a3", all = TRUE)
b <- c(1,round(0.5*max(agroforestry$freq, na.rm = TRUE)),max(agroforestry$freq, na.rm = TRUE))

p2 <- agroforestry %>%
  ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "#ceece7", #shade is 75% lighter using https://mdigi.tools/lighten-color/#aa4499
                      high = "#44AA99", 
                      na.value = "white",
                      limits = c(1,max(b)), 
                      breaks = b, 
                      labels = format(b)) + 
  ggtitle("B. Agroforestry")

#AVOIDED FOREST CONVERSION
avoidedFC <- fig1b %>% filter(pathway == "AFC")
avoidedFC <- merge(map, avoidedFC, by = "iso_a3", all = TRUE)
b <- c(1,round(0.5*max(avoidedFC$freq, na.rm = TRUE)),max(avoidedFC$freq, na.rm = TRUE))

p3 <- avoidedFC %>%
  ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "#eccee7", 
                      high = "#AA4499",
                      na.value = "white", 
                      limits = c(1,max(b)), 
                      breaks = b, 
                      labels = format(b)) + 
  ggtitle("C. Avoided Forest Conversion")

#COASTAL WETLAND
coastal <- fig1b %>% filter(pathway == "AWC/CWR")
coastal <- merge(map, coastal, by = "iso_a3", all = TRUE)
b <- c(1,round(0.5*max(coastal$freq, na.rm = TRUE)),max(coastal$freq, na.rm = TRUE))

p4 <- coastal %>%
  ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "#bcb0f4", 
                      high = "#392599", 
                      na.value = "white", 
                      limits = c(1,max(b)), 
                      breaks = b, 
                      labels = format(b)) + 
  ggtitle("D. Avoided Wetland Conversion & Coastal Wetland Restoration")

#GRASSLAND
grassland <- fig1b %>% filter(pathway == "AGC/GrR")
grassland <- merge(map, grassland, by = "iso_a3", all = TRUE)
b <- c(1,round(0.5*max(grassland$freq, na.rm = TRUE)),max(grassland$freq, na.rm = TRUE))

p5 <- grassland %>%
  ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "#fbf1c0", 
                      high = "#E6C20E", 
                      na.value = "white", 
                      limits = c(1,max(b)), breaks = b, labels = format(b)) + 
  ggtitle("E. Avoided Grassland Conversion & Grassland Restoration")

#PEATLAND
peatland <- fig1b %>% filter(pathway == "APC/PeR")
peatland <- merge(map, peatland, by = "iso_a3", all = TRUE)
b <- c(1,round(0.5*max(peatland$freq, na.rm = TRUE)),max(peatland$freq, na.rm = TRUE))

p6 <- peatland %>%
  ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "#ebd1f9", 
                      high = "#AF47E8", 
                      na.value = "white", 
                      limits = c(1,max(b)), 
                      breaks = b, 
                      labels = format(b)) +
  ggtitle("F. Avoided Peatland Conversion & Peatland Restoration")

#Reforestation (#fb8072)
reforestation <- fig1b %>% filter(pathway == "RFo")
reforestation <- merge(map, reforestation, by = "iso_a3", all = TRUE)
b <- c(1,round(0.5*max(reforestation$freq, na.rm = TRUE)),max(reforestation$freq, na.rm = TRUE))

p7 <- reforestation %>%
  ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "#fbddc0", 
                      high = "#E8790D", 
                      na.value = "white", 
                      limits = c(1,max(b)), 
                      breaks = b, labels = format(b)) + 
  ggtitle("G. Reforestation")

#CLIMATE SMART FORESTRY
climateSmart <- fig1b %>% filter(pathway == "CSF")
climateSmart <- merge(map, climateSmart, by = "iso_a3", all = TRUE)
b <- c(1,round(0.5*max(climateSmart$freq, na.rm = TRUE)),max(climateSmart$freq, na.rm = TRUE))

p8 <- climateSmart %>%
  ggplot() + 
  geom_sf(aes(fill = freq), position = "identity") + 
  labs(fill = 'Constraint Observations') +
  theme_map() + 
  theme(legend.title = element_text(size = 8), plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low = "#a5f8c0", 
                      high = "#09712B", 
                      na.value = "white", 
                      limits = c(1,max(b)), 
                      breaks = b, 
                      labels = format(b)) + 
  ggtitle("H. Climate Smart Forestry")

library(gridExtra)
jpeg("FIGURE1_2025.jpg", width = 1000, height = 1000)
grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, nrow = 4, ncol = 2)
dev.off()

```


##Figure 2A: Bar chart
```{r}
levels1 <- as.data.frame(table(df$constraintabb))
colnames(levels1) = c("constraintabb", "freq")

constraintCodes2 <- constraintCodes[,c("constraintabb", "constraintcategory")]
constraintCodes2 <- head(constraintCodes2, 39)


constraintCodes2 <- constraintCodes2 %>% pivot_wider(names_from = constraintabb, values_from = constraintcategory)
constraintCodes2 <- constraintCodes2 %>% pivot_longer(cols = everything(),values_to = "constraintcategory", names_to = "constraintabb") 

 constraintCodes2 <- constraintCodes2 %>%
   mutate(constraintabb = ifelse(constraintabb %in% c("blank1", "blank2", "blank3", "blank4", "blank5", "blank6", "blank7", "blank8"), NA, constraintabb))

levels1 <- right_join(levels1, constraintCodes2, by = "constraintabb")

levels1 <- levels1 %>% mutate(freq = ifelse(is.na(constraintabb) == TRUE, 0, freq))

levels <- levels1[with(levels1, order(levels1$constraintcategory, levels1$freq)),]

testLabels <- levels$constraintabb
testLabels <- testLabels[-c(47,48)]


levels1$constraintabb <- factor(levels1$constraintabb, levels$constraintabb) 



fig2 <- levels1  %>% filter(constraintcategory != "NA") %>% 
  filter(constraintabb != "NA") %>%
  group_by(constraintabb) %>% 
  mutate(constraintabb = fct_rev(constraintabb))%>%
  ggplot(aes(x = constraintabb, y = freq, fill = constraintcategory)) + 
  geom_bar(stat = 'identity') +
  scale_fill_manual(values= rev(c(

                             "#a65628",
                             "#e41a1c", 
                             "#ff7f00",
                             "#377eb8", 
                             "#984ea3",
                             "#ffff33",
                             "#4daf4a"))) + 
   scale_x_discrete(labels = rev(testLabels)) +
  coord_flip() + 
  theme_bw() + 
  theme(legend.title = element_text(size = 8)) +
  xlab("Constraint") + 
  ylab("Observations")   +
  labs(fill = "Constraint Category") + 
   geom_text(aes(label = freq), hjust = 1.25, colour = "lightgrey", size = 2.5)

#jpeg("FIGURE2.jpg", width = 800, height = 800)
fig2
#dev.off()


```



##Figure 2B: Pie chart of constraint categories
```{r}

constraint_counts <- df %>%
  group_by(constraintcategory) %>%
  summarise(total_count = n())

constraint_counts <- constraint_counts %>%
  mutate(percent = (total_count / sum(total_count)) * 100)
print(constraint_counts)


colors <- c("#4daf4a", "#ffff33", "#984ea3", "#377eb8", "#ff7f00",
            "#e41a1c", "#a65628") 


ggplot(constraint_counts, aes(x = "", y = total_count, fill = constraintcategory)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  scale_fill_manual(values = colors) +
  theme_void() +
  labs(fill = "Constraint Category", title = "Number of constraints in each category") +
  geom_text(aes(label = total_count), position = position_stack(vjust = 0.5), size = 6) +
  theme(legend.text = element_text(size = 12), legend.title = element_text(size = 14), plot.title = element_text(size = 16, face = "bold"))


```



##Figure 2C: Network of co-occuring constraints

```{r}

constraintabb_table <- df %>%
  distinct(constraintabb, Subregion) %>% 
  mutate(presence = 1) %>% 
  spread(key = Subregion, value = presence, fill = 0) 

constraintabb_binary <- constraintabb_table[,-1]

jaccard_distance_matrix <- vegdist(constraintabb_binary, method = "jaccard")
jaccard_similarity_matrix <- 1 - jaccard_distance_matrix
jaccard_similarity_df <- as.data.frame(as.matrix(jaccard_similarity_matrix))

rownames(jaccard_similarity_df) <- constraintabb_table$constraintabb
colnames(jaccard_similarity_df) <- constraintabb_table$constraintabb

jaccard_similarity_long <- jaccard_similarity_df %>%
  mutate(Constraint1 = rownames(.)) %>%
  gather(key = "Constraint2", value = "Similarity", -Constraint1)

jaccard_similarity_long <- jaccard_similarity_long %>%
  filter(Constraint1 != Constraint2)

average_jaccard_similarity <- mean(jaccard_similarity_long$Similarity)
print(paste("Average Jaccard Similarity:", average_jaccard_similarity))


highest_similar_pairs <- lapply(1:nrow(jaccard_similarity_df), function(i) {
  max_val <- max(jaccard_similarity_df[i, -i]) # Exclude self-similarity (1)
  if (max_val > 0.1) {
    max_index <- which(jaccard_similarity_df[i, ] == max_val)
    return(c(i, max_index))
  } else {
    return(NA)
  }
})

highest_similar_pairs <- Filter(Negate(is.na), highest_similar_pairs)

constraint1_indices <- sapply(highest_similar_pairs, `[`, 1)
constraint2_indices <- sapply(highest_similar_pairs, `[`, 2)

# Dataframe with the most similar pairs and their similarity values
similar_pairs_df <- data.frame(
  Constraint1 = rownames(jaccard_similarity_df)[unlist(constraint1_indices)],
  Constraint2 = rownames(jaccard_similarity_df)[unlist(constraint2_indices)],
  Similarity = mapply(function(x, y) jaccard_similarity_df[x, y], unlist(constraint1_indices), unlist(constraint2_indices))
)

similar_pairs_df <- similar_pairs_df %>%
  rowwise() %>%
  mutate(
    SortedPair = paste(sort(c(Constraint1, Constraint2)), collapse = "-")
  ) %>%
  ungroup() %>%
  filter(!duplicated(SortedPair)) %>%
  select(-SortedPair)

similar_pairs_with_buckets <- similar_pairs_df %>%
  left_join(constraintCodes, by = c("Constraint1" = "constraintabb")) %>%
  rename(Constraint1Bucket = constraintcategory) %>%
  left_join(constraintCodes, by = c("Constraint2" = "constraintabb")) %>%
  rename(Constraint2Bucket = constraintcategory)

same_bucket_count <- nrow(similar_pairs_with_buckets %>%
                            filter(Constraint1Bucket == Constraint2Bucket))

different_bucket_count <- nrow(similar_pairs_with_buckets %>%
                                 filter(Constraint1Bucket != Constraint2Bucket))

# Print the results
print(paste("Number of pairs from the same category:", same_bucket_count))
print(paste("Number of pairs from different category:", different_bucket_count))


edges_df <- data.frame(
  from = similar_pairs_df$Constraint1,
  to = similar_pairs_df$Constraint2,
  weight = similar_pairs_df$Similarity
)

edges_df <- edges_df[edges_df$weight > 0.1, ]  # Apply a threshold to filter weak connections


g <- graph_from_data_frame(edges_df, directed = FALSE)

vertex_data <- data.frame(name = V(g)$name)
vertex_data <- merge(vertex_data, constraintCodes, by.x = "name", by.y = "constraintabb", all.x = TRUE)
vertex_data$constraintcategory <- factor(vertex_data$constraintcategory, levels = unique(constraintCodes$constraintcategory))


color_mapping <- c(
  "Economic" = "#4daf4a",
  "Governments & Organizations" = "#ffff33",
  "Knowledge" = "#984ea3",
  "Material Inputs" = "#377eb8",
  "Policies & Rules" = "#e41a1c",
  "Social & Behavioral" = "#a65628",
  "Negative Side Effects" = "#ff7f00"
)

g_tidy <- as_tbl_graph(g) %>%
  left_join(vertex_data, by = c("name" = "name"))

ggraph(g_tidy, layout = "fr") + # Using Fruchterman-Reingold layout
  geom_edge_link(aes(edge_alpha = weight, edge_width = weight), color = "gray4", show.legend = FALSE) +
  geom_node_point(aes(color = constraintcategory), size = 8) +
  geom_node_text(aes(label = name), repel = TRUE) +  
  scale_color_manual(values = color_mapping) +
  theme_void() +
  labs(title = "Constraint with Highest Co-occurrence within Subregion") +
  theme(legend.position = "right") +
  scale_edge_width(range = c(0.5, 2))  

```



##Figure 3: Constraint by pathway
```{r}
fig3 <- as.data.frame(table(df$constraintabb, df$pathway))
colnames(fig3) <- c("constraintabb", "pathway", "freq")

p1 <- fig3  %>% filter(pathway != "NA") %>% 
  group_by(constraintabb) %>% 
   mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = constraintabb, y = freq, fill = pathway)) + 
  geom_col(position =  position_fill(reverse = TRUE)) +
  scale_fill_manual(values = pathwayColors) + 
  coord_flip() + 
  theme_bw() + 
  theme(legend.position = "none") +
  xlab("Constraint") +
    ylab("Proportion of Constraint Observations")
  
p2 <- fig3 %>% 
  filter(pathway != "NA") %>% 
  group_by(pathway) %>% 
  mutate(freq2 = freq/sum(freq)) %>% 
  mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = pathway, y = constraintabb, fill = pathway)) + 
  geom_tile(aes(alpha = freq2), color = "white") + 
  theme(axis.text.x = element_text(angle = 90)) + 
  scale_fill_manual(values = pathwayColors) +
  theme_classic()+
  theme(axis.text.x=element_text(angle=0,hjust=0.5,vjust=0.5),
        axis.text.y=element_blank(), axis.ticks.y = element_blank(),
        axis.title.y= element_blank()) +
  geom_text(aes(label = paste0(round(freq2,2)*100,"%")), size = 3, nudge_x = -0.1) + 
  geom_text(aes(label = paste0("(",round(freq,2),")")),size = 2, nudge_x = 0.1, color = "#5A5A5A") +
  scale_alpha(guide = 'none') + 
  xlab("Pathway")

grid.arrange(p1, p2, ncol = 2)

```


##Figure 4: Constraint by SDG region
```{r}
fig4 <- as.data.frame(table(df$constraintabb, df$SDGrgion))
colnames(fig4) <- c("constraintabb", "region", "freq")

p1 <- fig4  %>% filter(region != "NA") %>% 
  group_by(constraintabb) %>% 
  mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = constraintabb, y = freq, fill = region)) + 
  geom_col(position = position_fill(reverse = TRUE)) +
  scale_fill_manual(values = sdgColors) + 
  coord_flip() + 
  theme_bw() + 
  theme(legend.position = "none") +
  xlab("Constraint") + 
  ylab("Proportion of Constraint Observations")

p2 <- fig4 %>% 
  filter(region != "NA") %>% 
  group_by(region) %>% 
  mutate(freq2 = freq/sum(freq)) %>% 
  mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = region, y = constraintabb, fill = region)) + 
  geom_tile(aes(alpha = freq2), color = "white") + 
  theme(axis.text.x = element_text(angle = 90)) + 
  scale_fill_manual(values = sdgColors) +
  theme_classic()+
  theme(axis.text.y=element_blank(), 
        axis.ticks.y = element_blank(),
        axis.title.y= element_blank()) +
  scale_x_discrete(labels = c("CSA", "ESEA", "ENA", "LAC", "NAWA", "OCA","SSA")) +
  geom_text(aes(label= round(freq, 2)), size = 2) + 
  scale_alpha(guide = 'none') + 
  xlab("SDG Region")

p2 <- fig4 %>% 
  filter(region != "NA") %>% 
  group_by(region) %>% 
  mutate(freq2 = freq/sum(freq)) %>% 
  mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = region, y = constraintabb, fill = region)) + 
  geom_tile(aes(alpha = freq2), color = "white") + 
  theme(axis.text.x = element_text(angle = 90)) + 
  scale_fill_manual(values = sdgColors) +
  theme_classic()+
  theme(axis.text.x=element_text(angle=0,hjust=0.5,vjust=0.5),
        axis.text.y=element_blank(), axis.ticks.y = element_blank(),
        axis.title.y= element_blank()) +
  scale_x_discrete(labels = c("CSA", "ESEA", "ENA", "LAC", "NAWA", "OCA","SSA")) +
  geom_text(aes(label = paste0(round(freq2,2)*100,"%")), size = 3, nudge_x = -0.1) + 
  geom_text(aes(label = paste0("(",round(freq,2),")")),size = 2, nudge_x = 0.1, color = "#5A5A5A") +
  scale_alpha(guide = 'none') + 
  xlab("SDG Region")

grid.arrange(p1, p2, ncol = 2)

```

##Figure 5: Maps of top constraint in UN subregions
```{r}

unique_pathways <- unique(df$pathway)

pathway_dataframes <- list()

for (path in unique_pathways) {
  assign(paste("df_", path, sep = ""), df[df$pathway == path, ])
}

df_AGC <- `df_AGC/GrR`
df_APC <- `df_APC/PeR`
df_AWC <- `df_AWC/CWR`

dataframe_names <- c("df_AFC", "df_AGC", "df_AgFo", "df_APC", "df_AWC", "df_CSF", "df_RFo", "df")


process_data <- function(input_df) {
  df_summary <- input_df %>%
    group_by(Subregion) %>%
    count(constraintabb) %>%
    top_n(1, n) %>%
    ungroup()
  
  df_summary_updated <- df_summary %>%
    add_count(Subregion, name = "count") %>%
    group_by(Subregion) %>%
    mutate(
      constraintabb = ifelse(count > 1, "Tie", constraintabb),
      n = ifelse(count > 1, NA, n)
    ) %>%
    filter(row_number() == 1) %>%
    select(-count)
  
  world <- ne_countries(scale = "medium", returnclass = "sf")
  world_data <- left_join(world, df_summary_updated, by = c("subregion" = "Subregion"))
  
  return(world_data)
}


df_RFo_processed <- process_data(df_RFo)
df_AgFo_processed <- process_data(df_AgFo)
df_AFC_processed <- process_data(df_AFC)
df_processed <- process_data(df)


df_RFo_processed <- df_RFo_processed %>%
  mutate(constraintabb = str_trim(constraintabb))

df_AgFo_processed <- df_AgFo_processed %>%
  mutate(constraintabb = str_trim(constraintabb))


df_AFC_processed <- df_AFC_processed %>%
  mutate(constraintabb = str_trim(constraintabb))

df_processed <- df_processed %>%
  mutate(constraintabb = str_trim(constraintabb))


names(constraint_colors) <- str_trim(names(constraint_colors))


plot_all_pathways <- ggplot(data = df_processed) +
  geom_sf(aes(fill = constraintabb), color = "white", size = 0.2) +
  scale_fill_manual(values = constraint_colors, na.value = "grey90") + 
  labs(title = "Most Common Constraints by SDG Subregion",
       fill = "Constraint") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

plot_Reforestation <- ggplot(data = df_RFo_processed) +
  geom_sf(aes(fill = constraintabb), color = "white", size = 0.2) +
  scale_fill_manual(values = constraint_colors, na.value = "grey90") +  
  labs(title = "Most Common Constraints by SDG Subregion: Reforestation",
       fill = "Constraint") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


plot_Agroforestry <- ggplot(data = df_AgFo_processed) +
  geom_sf(aes(fill = constraintabb), color = "white", size = 0.2) +
  scale_fill_manual(values = constraint_colors, na.value = "grey90") +  
  labs(title = "Most Common Constraints by SDG Subregion: Agroforestry",
       fill = "Constraint") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

plot_AFC <- ggplot(data = df_AFC_processed) +
  geom_sf(aes(fill = constraintabb), color = "white", size = 0.2) +
  scale_fill_manual(values = constraint_colors, na.value = "grey90") +  
  labs(title = "Most Common Constraints by SDG Subregion: Avoided Forest Conversion",
       fill = "Constraint") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())




plot_Reforestation
plot_Agroforestry
plot_AFC
plot_all_pathways

```

##Figure 5: Pie charts for tied subregions

```{r}
generate_pie_charts <- function(df, title_suffix) {
  ties <- df %>%
    group_by(Subregion, constraintabb) %>%
    summarize(count = n(), .groups = "drop") %>%
    group_by(Subregion) %>%
    mutate(max_count = max(count)) %>%
    filter(count == max_count) %>%
    ungroup()

  subregion_with_ties <- ties %>%
    add_count(Subregion, name = "num_ties") %>%
    filter(num_ties > 1) %>%
    select(Subregion) %>%
    distinct()

  top_tied_constraints <- ties %>%
    filter(Subregion %in% subregion_with_ties$Subregion)

  top_tied_constraints <- top_tied_constraints %>%
    group_by(Subregion) %>%
    mutate(perc = count / sum(count)) %>%
    ungroup()

  ggplot(top_tied_constraints, aes(x="", y=perc, fill=constraintabb)) +
    geom_bar(stat="identity", width=1) +
    coord_polar("y", start=0) +
    facet_wrap(~Subregion) +
    theme_void() +
    scale_fill_manual(values = constraint_colors) +
    theme(legend.title = element_blank(),
          strip.text = element_text(size = 16, face = "bold")) +  
    labs(title = paste("Top Tied Constraints per Subregion for", title_suffix),
        y = "Percentage") +
    guides(fill=guide_legend(title="Constraintabb")) +
    scale_y_continuous(labels = scales::percent)

}

#Pie charts for each pathway 

management_pie_charts <- generate_pie_charts(df_AgFo, "Agroforestry")
management_pie_charts

restoration_pie_charts <- generate_pie_charts(df_RFo, "Reforestation")
restoration_pie_charts

protection_pie_charts <- generate_pie_charts(df_AFC, "Avoided Forest Converstion")
protection_pie_charts

all_pie_charts <- generate_pie_charts(df, "All Pathways")
all_pie_charts
```




## Figure 5: Legend

```{r}
#First find the constraints to include in the legend (those included in the maps and/or pie charts)


generate_pie_charts <- function(df, title_suffix) {
  # Identify subregions with ties and the corresponding top constraints
  ties <- df %>%
    group_by(Subregion, constraintabb) %>%
    summarize(count = n(), .groups = "drop") %>%
    group_by(Subregion) %>%
    mutate(max_count = max(count)) %>%
    filter(count == max_count) %>%
    ungroup()

  subregion_with_ties <- ties %>%
    add_count(Subregion, name = "num_ties") %>%
    filter(num_ties > 1) %>%
    select(Subregion) %>%
    distinct()

  # Get all constraints that are tied for the top spot 
  top_tied_constraints <- ties %>%
    filter(Subregion %in% subregion_with_ties$Subregion)

  return(unique(top_tied_constraints$constraintabb))
}

unique_management_constraints <- generate_pie_charts(df_AgFo, "Agroforestry")
unique_restoration_constraints <- generate_pie_charts(df_RFo, "Reforestation")
unique_protection_constraints <- generate_pie_charts(df_AFC, "Avoided Forest Conversion")
unique_all_constraints <- generate_pie_charts(df, "All Pathways")

unique_constraints_pie_charts <- unique(c(
  unique_management_constraints, 
  unique_restoration_constraints, 
  unique_protection_constraints, 
  unique_all_constraints
))

unique_RFo_constraints <- unique(df_RFo_processed$constraintabb)
unique_AgFo_constraints <- unique(df_AgFo_processed$constraintabb)
unique_AFC_constraints <- unique(df_AFC_processed$constraintabb)
unique_all_constraints_maps <- unique(df_processed$constraintabb)

unique_constraints_maps <- unique(c(
  unique_RFo_constraints, 
  unique_AgFo_constraints, 
  unique_AFC_constraints, 
  unique_all_constraints_maps
))

all_unique_constraints <- unique(c(unique_constraints_pie_charts, unique_constraints_maps))

all_unique_constraints <- all_unique_constraints[!is.na(all_unique_constraints)]



all_constraints <- c(
  "Lack of funding",
  "Lack of credit",
  "Lack of insurance",
  "Inadequate markets for ecosystem services",
  "Inadequate markets for outputs",
  "Inadequate prices for ecosystem services",
  "Inadequate prices for outputs",
  "Time lag in benefits",
  "Tradeoffs with agriculture",
  "Insufficient info about yields or profits",
  "Lack of info about co-benefits",
  "Lack of info on how to design or manage",
  "Lack of or insufficient technical support",
  "Limited land manager capability",
  "Inadequate planning or management",
  "Preferences for non-NCS land uses",
  "Disinterest or skepticism",
  "Limited social learning/exchange networks",
  "Lack of opportunity to participate",
  "Interpersonal conflict",
  "Human-wildlife conflict",
  "Concerns over negative equity impacts",
  "Ineffective laws, policies, or regulations",
  "Lack of laws, policies, or regulations",
  "Insecure land tenure",
  "Unsecure or uncertain benefits",
  "Incentives for non-NCS",
  "Lack of administrative capacity",
  "Lack of enforcement",
  "Lack of coordination among organizations",
  "Lack of political will",
  "Violent conflict or civil unrest",
  "Corruption or lack of transparency",
  "Lack of input materials or infrastructure",
  "Lack of labor",
  "Lack of water/water distribution networks",
  "Unsuitable land",
  "Tradeoffs with other land uses",
  "Other negative side effects",
  "NA"
)

constraints_present <- all_constraints[all_constraints %in% all_unique_constraints]

constraints_absent <- setdiff(all_constraints, constraints_present)

cat("Constraints that appear in figures:\n")
print(constraints_present)

cat("\nConstraints that do NOT appear in any figures:\n")
print(constraints_absent)



constraint_colors_filtered <- constraint_colors[names(constraint_colors) %in% constraints_present]

constraintCodes2_filtered <- constraintCodes2[constraintCodes2$constraintabb %in% constraints_present, ]

constraints_df <- data.frame(
  Constraint = names(constraint_colors_filtered),
  Color = constraint_colors_filtered
)

constraints_df <- merge(constraints_df, constraintCodes2_filtered, 
                        by.x = "Constraint", by.y = "constraintabb", all.x = TRUE)

default_color <- "#D3D3D3" 
constraints_df$Color[is.na(constraints_df$Color)] <- default_color

ordered_categories <- sort(unique(constraints_df$constraintcategory[!constraints_df$constraintcategory %in% c("Tie", "NA")]))
constraints_df$constraintcategory <- factor(
  constraints_df$constraintcategory, 
  levels = c(ordered_categories, "Tie", "NA")
)

constraints_df <- constraints_df[order(constraints_df$constraintcategory, constraints_df$Constraint), ]

create_legend_grob <- function(df) {
  grobs <- lapply(1:nrow(df), function(i) {
    legend_color <- rectGrob(width = unit(1, "lines"), height = unit(1, "lines"), 
                             gp = gpar(fill = df$Color[i], col = "black"))
    legend_text <- textGrob(df$Constraint[i], x = unit(1.5, "lines"), just = "left")
    legend_item <- gTree(children = gList(legend_color, legend_text), cl = "legend_item")
    return(legend_item)
  })
  return(grobs)
}

legend_grobs <- list()
for (category in levels(constraints_df$constraintcategory)) {
  category_df <- subset(constraints_df, constraintcategory == category)
  
  if (nrow(category_df) > 0) { 
    category_title <- textGrob(category, gp = gpar(fontface = "bold"), just = "left", x = unit(0, "npc"))
    
    category_grobs <- create_legend_grob(category_df)
    
    legend_grobs <- c(legend_grobs, list(category_title), category_grobs, list(textGrob(" ")))  # Add space between categories
  }
}

if (length(legend_grobs) > 0) {
  legend_grobs <- legend_grobs[-length(legend_grobs)]
}

legend_grob <- arrangeGrob(
  grobs = do.call(gList, legend_grobs), 
  ncol = 1, 
  vp = viewport(x = 0.5, y = 0.5, just = c("center", "center"))
)

grid.newpage()
grid.draw(legend_grob)

```


##Figure 6A: Constraints in 5 countries with highest Climate Mitigation Potential
```{r}

Naturebase_sorted <- Naturebase[order(Naturebase$total_ncs_co2e_total, decreasing = TRUE), ]
top_5_countries <- head(Naturebase_sorted, 5)$ISO
top_5_countries

top_countries <- df[df$ISOCode %in% top_5_countries, ]

top_countries <- merge(top_countries, Naturebase[c("ISOCode", "total_ncs_co2e_total")], by = "ISOCode")
names(top_countries)[names(top_countries) == "total_ncs_co2e_total"] <- "carbon_potential"

top_countries_clean <- top_countries


country_totals <- Naturebase %>%
                  group_by(ISOCode) %>%
                  summarize(Total = sum(total_ncs_co2e_total, na.rm = TRUE))

top_countries <- merge(top_countries, country_totals, by = "ISOCode")

top_countries$ISOCode <- with(top_countries, reorder(ISOCode, -Total))


bucket_colors <- rev(c("#4daf4a", "#ffff33","#984ea3",  "#377eb8","#ff7f00", "#e41a1c",  "#a65628"))

top_countries$ISOCode <- with(top_countries, reorder(ISOCode, Total))

total_constraints_per_country <- top_countries %>%
  group_by(ISOCode) %>%
  summarize(TotalConstraints = n())

top_countries <- merge(top_countries, total_constraints_per_country, by = "ISOCode")

top_countries <- top_countries %>%
  group_by(ISOCode, constraintcategory) %>%
  summarize(Count = n(), .groups = 'drop')

top_countries <- merge(top_countries, total_constraints_per_country, by = "ISOCode")

top_countries <- merge(top_countries, placeCodes, by = "ISOCode")

top_countries <- merge(top_countries, country_totals, by = "ISOCode")

top_countries <- top_countries %>%
  mutate(Percentage = (Count / TotalConstraints)*100)

top_countries$Country <- with(top_countries, paste(Country, "(n =", TotalConstraints, ")"))

top_countries$Country <- with(top_countries, reorder(Country, Total))



top_countries$constraintcategory <- factor(top_countries$constraintcategory, levels = rev(levels(factor(top_countries$constraintcategory))))

fig6A <- ggplot(top_countries, aes(x = Country, y = Percentage, fill = constraintcategory)) +
  geom_bar(stat = "identity", position = "fill") +
  scale_fill_manual(values = bucket_colors) +  
  labs(title = "Proportion of Constraints by Country and Constraint Bucket",
       x = "Country", 
       y = "Proportion of Constraints",
       fill = "Constraint Types") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        legend.title = element_text(size = 8)) +
  coord_flip()

fig6A

```



##Figure 6B: Top constraint by country map

```{r}

df <- df %>%
  mutate(constraintabb = str_trim(constraintabb))


process_data <- function(input_df) {
  df_summary <- input_df %>%
    group_by(ISOCode) %>%
    count(constraintabb) %>%
    top_n(1, n) %>%
    ungroup()
  
  df_summary_updated <- df_summary %>%
    add_count(ISOCode, name = "count") %>%
    group_by(ISOCode) %>%
    mutate(
      constraintabb = ifelse(count > 1, "Tie", constraintabb),
      n = ifelse(count > 1, NA, n)
    ) %>%
    filter(row_number() == 1) %>%
    select(-count)
  
  world <- ne_countries(scale = "medium", returnclass = "sf")
  world_data <- left_join(world, df_summary_updated, by = c("iso_a3" = "ISOCode"))
  
  return(world_data)
}


topcountraint_country <- process_data(df)


p_6B <- ggplot(data = topcountraint_country) +
  geom_sf(aes(fill = constraintabb), color = "white", size = 0.2) +
  scale_fill_manual(values = constraint_colors, na.value = "grey90") +  
  labs(title = "Most Common Constraints by Country",
       fill = "Constraint") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

p_6B

```

##Figure 6B: Legend

```{r}
#Identify the constraints included in Figure 6B and integrate them in a legend

constraints_in_map <- unique(topcountraint_country$constraintabb)

constraints_in_map <- constraints_in_map[!is.na(constraints_in_map)]

constraint_colors_filtered <- constraint_colors[names(constraint_colors) %in% constraints_in_map]

constraintCodes2_filtered <- constraintCodes2[constraintCodes2$constraintabb %in% constraints_in_map, ]

constraints_df <- data.frame(
  Constraint = names(constraint_colors_filtered),
  Color = constraint_colors_filtered
)

constraints_df <- merge(constraints_df, constraintCodes2_filtered, 
                        by.x = "Constraint", by.y = "constraintabb", all.x = TRUE)

default_color <- "#D3D3D3"  
constraints_df$Color[is.na(constraints_df$Color)] <- default_color

ordered_categories <- sort(unique(constraints_df$constraintcategory[!constraints_df$constraintcategory %in% c("Tie", "NA")]))
constraints_df$constraintcategory <- factor(
  constraints_df$constraintcategory, 
  levels = c(ordered_categories, "Tie", "NA")
)

constraints_df <- constraints_df[order(constraints_df$constraintcategory, constraints_df$Constraint), ]

create_legend_grob <- function(df) {
  grobs <- lapply(1:nrow(df), function(i) {
    legend_color <- rectGrob(width = unit(1, "lines"), height = unit(1, "lines"), 
                             gp = gpar(fill = df$Color[i], col = "black"))
    legend_text <- textGrob(df$Constraint[i], x = unit(1.5, "lines"), just = "left")
    legend_item <- gTree(children = gList(legend_color, legend_text), cl = "legend_item")
    return(legend_item)
  })
  return(grobs)
}

legend_grobs <- list()
for (category in levels(constraints_df$constraintcategory)) {
  category_df <- subset(constraints_df, constraintcategory == category)
  
  if (nrow(category_df) > 0) {  
    category_title <- textGrob(category, gp = gpar(fontface = "bold"), just = "left", x = unit(0, "npc"))

    category_grobs <- create_legend_grob(category_df)

    legend_grobs <- c(legend_grobs, list(category_title), category_grobs, list(textGrob(" ")))  # Add space between categories
  }
}


if (length(legend_grobs) > 0) {
  legend_grobs <- legend_grobs[-length(legend_grobs)]
}


legend_grob <- arrangeGrob(
  grobs = do.call(gList, legend_grobs), 
  ncol = 1, 
  vp = viewport(x = 0.5, y = 0.5, just = c("center", "center"))
)


grid.newpage()
grid.draw(legend_grob)


```


#Supplemental Information Figures


##Figure S2: Constraints by Income
```{r}
figS6 <- as.data.frame(table(df$constraintabb, df$Income))
colnames(figS6) <- c("constraintabb", "income", "freq")

figS6$income <- factor(figS6$income, levels = c("low", "lower-middle", "upper-middle", "high"))


p1 <- figS6 %>% 
  filter(income != "NA-Income") %>% 
  group_by(constraintabb) %>% 
  mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = constraintabb, y = freq, fill = income)) + 
  geom_col(position = position_fill(reverse = TRUE)) +
  scale_fill_manual(values = incomeColors) + 
  coord_flip() + 
  theme_bw() + 
  theme(legend.position = "none") +
  xlab("Constraint") + 
  ylab("Proportion of Constraint Observations")


p2 <- figS6 %>% 
  filter(income != "NA-Income") %>% 
  group_by(income) %>% 
  mutate(freq2 = freq/sum(freq)) %>% 
  mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = income, y = constraintabb, fill = income)) + 
  geom_tile(aes(alpha = freq2), color = "white") + 
  theme(axis.text.x = element_text(angle = 90)) + 
  scale_fill_manual(values = incomeColors) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5),
        axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        axis.title.y = element_blank()) +
  geom_text(aes(label = paste0(round(freq2, 2) * 100, "%")), size = 3, nudge_x = -0.1) + 
  geom_text(aes(label = paste0("(", round(freq, 2), ")")), size = 2, nudge_x = 0.1, color = "#5A5A5A") +
  scale_alpha(guide = 'none') + 
  xlab("Income")

# Arrange plots
grid.arrange(p1, p2, ncol = 2)

```


##Figure S3: Constraint ranking by countries
```{r}
levels1 <- df %>%
  group_by(constraintabb) %>%
  summarize(freq = n_distinct(ISOCode), .groups = "drop")

colnames(levels1) <- c("constraintabb", "freq")

constraintCodes2 <- constraintCodes[, c("constraintabb", "constraintcategory")]
constraintCodes2 <- head(constraintCodes2, 39)

constraintCodes2 <- constraintCodes2 %>%
  pivot_wider(names_from = constraintabb, values_from = constraintcategory) %>%
  pivot_longer(cols = everything(), values_to = "constraintcategory", names_to = "constraintabb")

constraintCodes2 <- constraintCodes2 %>%
  mutate(constraintabb = ifelse(constraintabb %in% c("blank1", "blank2", "blank3", "blank4", "blank5", "blank6", "blank7", "blank8"), NA, constraintabb))

levels1 <- right_join(levels1, constraintCodes2, by = "constraintabb")

levels1 <- levels1 %>%
  mutate(freq = ifelse(is.na(constraintabb), 0, freq))

levels <- levels1 %>%
  arrange(constraintcategory, freq)

testLabels <- levels$constraintabb
testLabels <- testLabels[-c(47, 48)]

levels1$constraintabb <- factor(levels1$constraintabb, levels$constraintabb)

figs3 <- levels1 %>%
  filter(!is.na(constraintcategory)) %>%
  filter(!is.na(constraintabb)) %>%
  group_by(constraintabb) %>%
  mutate(constraintabb = fct_rev(constraintabb)) %>%
  ggplot(aes(x = constraintabb, y = freq, fill = constraintcategory)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = rev(c(
                             "#a65628",
                             "#e41a1c", 
                             "#ff7f00",
                             "#377eb8", 
                             "#984ea3",
                             "#ffff33",
                             "#4daf4a"
  ))) +
  scale_x_discrete(labels = rev(testLabels)) +
  coord_flip() +
  theme_bw() +
  theme(legend.title = element_text(size = 8)) +
  xlab("Constraint") +
  ylab("Unique Countries Observed") +
  labs(fill = "Constraint Category") +
  geom_text(aes(label = freq), hjust = 1.25, colour = "lightgrey", size = 2.5)

figs3

```


##Figure S4: Constraints rankings in countries with top 10% highest CMP
```{r}

total_co2e_sum <- sum(Naturebase$total_ncs_co2e_total)

Naturebase$percent_total_co2e <- (Naturebase$total_ncs_co2e_total / total_co2e_sum) * 100

Naturebase$percent_total_co2e <- (Naturebase$total_ncs_co2e_total / total_co2e_sum) * 100

cutoff <- quantile(Naturebase$percent_total_co2e, 0.90)

higher_percent_countries <- Naturebase[Naturebase$percent_total_co2e >= cutoff, ]

higher_percent_sum <- sum(higher_percent_countries$percent_total_co2e)


filtered_data <- df[df$ISOCode %in% higher_percent_countries$ISOCode, ]

levels1 <- as.data.frame(table(filtered_data$constraintabb))
colnames(levels1) = c("constraintabb", "freq")

constraintCodes2 <- constraintCodes[,c("constraintabb", "constraintcategory")]
constraintCodes2 <- head(constraintCodes2, 39)

constraintCodes2 <- constraintCodes2 %>% pivot_wider(names_from = constraintabb, values_from = constraintcategory)
constraintCodes2 <- constraintCodes2 %>% pivot_longer(cols = everything(),values_to = "constraintcategory", names_to = "constraintabb")                                                

levels1 <- right_join(levels1, constraintCodes2, by = "constraintabb")

levels1 <- levels1 %>% mutate(freq = ifelse(is.na(constraintabb) == TRUE, 0, freq))

levels <- levels1[with(levels1, order(levels1$constraintcategory, levels1$freq)),]

testLabels <- levels$constraintabb
testLabels <- testLabels[-c(47,48)]


levels1$constraintabb <- factor(levels1$constraintabb, levels$constraintabb) 

figS4 <- levels1  %>% filter(constraintcategory != "NA") %>% 
  filter(constraintabb != "NA") %>%
  group_by(constraintabb) %>% 
  mutate(constraintabb = fct_rev(constraintabb))%>%
  ggplot(aes(x = constraintabb, y = freq, fill = constraintcategory)) + 
  geom_bar(stat = 'identity') +
  scale_fill_manual(values= rev(c(
                             "#a65628",
                             "#e41a1c",
                             "#ff7f00",
                             "#377eb8", 
                             "#984ea3",
                             "#ffff33",
                             "#4daf4a"))) + 
   scale_x_discrete(labels = rev(testLabels)) +
  coord_flip() + 
  theme_bw() + 
  theme(legend.title = element_text(size = 8)) +
  xlab("Constraint") + 
  ylab("Observations from Top 10% CSP Countries")   +
  labs(fill = "Constraint Category") + 
   geom_text(aes(label = freq), hjust = 1.25, colour = "lightgrey", size = 2.5)

figS4


```

##Figure S5: Constraint rankings for 5 countries with most data
```{r}

top5_countries <- df %>% 
  group_by(Country) %>% 
  summarise(n = n(), .groups = "drop") %>% 
  arrange(desc(n)) %>% 
  slice(1:5) %>% 
  pull(Country)

df <- df %>% 
  mutate(panel = ifelse(Country %in% top5_countries, Country, "Other"))

levels_panel <- df %>%
  group_by(panel, constraintabb) %>%
  summarise(freq = n(), .groups = "drop")

constraintCodes2 <- constraintCodes[, c("constraintabb", "constraintcategory")]
constraintCodes2 <- head(constraintCodes2, 39)
constraintCodes2 <- constraintCodes2 %>% 
  pivot_wider(names_from = constraintabb, values_from = constraintcategory) %>%
  pivot_longer(cols = everything(), values_to = "constraintcategory", names_to = "constraintabb") %>%
  mutate(constraintabb = ifelse(constraintabb %in% c("blank1", "blank2", "blank3", 
                                                     "blank4", "blank5", "blank6", 
                                                     "blank7", "blank8"), NA, constraintabb))

levels_panel <- right_join(levels_panel, constraintCodes2, by = "constraintabb") %>%
  mutate(freq = ifelse(is.na(freq), 0, freq))

levels_order <- levels_panel %>% 
  arrange(constraintcategory, freq) %>% 
  distinct(constraintabb, .keep_all = TRUE)
levels_panel$constraintabb <- factor(levels_panel$constraintabb, levels = rev(levels_order$constraintabb))

all_panels <- unique(df$panel)
all_constraints <- levels_order %>% select(constraintabb, constraintcategory)
levels_panel <- expand.grid(panel = all_panels, constraintabb = levels_order$constraintabb, 
                            stringsAsFactors = FALSE) %>%
  left_join(all_constraints, by = "constraintabb") %>%
  left_join(levels_panel %>% select(panel, constraintabb, freq), by = c("panel", "constraintabb"))
levels_panel$freq[is.na(levels_panel$freq)] <- 0
levels_panel$panel <- factor(levels_panel$panel, levels = c(top5_countries, "Other"))
levels_panel$constraintabb <- factor(levels_panel$constraintabb, levels = rev(levels_order$constraintabb))


desired_category_order <- c("Economic", "Governments & Organizations", "Knowledge", 
                            "Material Inputs", "Negative Side Effects", "Policies & Rules", 
                            "Social & Behavioral")
levels_panel <- levels_panel %>%
  mutate(constraintcategory = factor(constraintcategory, levels = desired_category_order))


plot_panel <- function(panel_name, show_y = TRUE) {
  data_panel <- levels_panel %>% 
    filter(panel == panel_name,
           !is.na(constraintabb),
           !is.na(constraintcategory))
  
  max_freq <- max(data_panel$freq, na.rm = TRUE)

  if(panel_name == "Other"){
    tick_step <- if(max_freq < 50) 1 else 50
  } else {
    tick_step <- if(max_freq < 5) 1 else 5
  }
  tick_breaks <- seq(0, max_freq, by = tick_step)
  
  p <- ggplot(data_panel, aes(x = constraintabb, y = freq, fill = constraintcategory)) +
    geom_bar(stat = "identity") +
    scale_fill_manual(values = c("#4daf4a", "#ffff33", "#984ea3", 
                             "#377eb8", "#ff7f00", "#e41a1c", "#a65628"))+
    coord_flip() +
    theme_bw() +
    xlab("Constraint") +
    ylab("Observations") +
    labs(fill = "Constraint Category") +
    geom_text(aes(label = freq), hjust = 1.25, colour = "lightgrey", size = 2.5) +
    ggtitle(panel_name) +
    scale_y_continuous(breaks = tick_breaks)
  
  if (!show_y) {
    p <- p + theme(axis.text.y = element_blank(),
                   axis.ticks.y = element_blank())
  }
  return(p)
}


p1 <- plot_panel(top5_countries[1], show_y = TRUE)    
p2 <- plot_panel(top5_countries[2], show_y = FALSE)
p3 <- plot_panel(top5_countries[3], show_y = FALSE)
p4 <- plot_panel(top5_countries[4], show_y = TRUE)    
p5 <- plot_panel(top5_countries[5], show_y = FALSE)
p6 <- plot_panel("Other", show_y = FALSE)

figS5 <- (p1 + p2 + p3) / (p4 + p5 + p6) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

figS5


```

##Figure S6: Constraint category breakdowns for 5 countries with most data
```{r}

top5_countries <- df %>% 
  group_by(Country) %>% 
  summarise(n = n(), .groups = "drop") %>% 
  arrange(desc(n)) %>% 
  slice(1:5) %>% 
  pull(Country)

df <- df %>% 
  mutate(panel = ifelse(Country %in% top5_countries, Country, "Other"))


levels_panel <- df %>%
  group_by(panel, constraintabb) %>%
  summarise(freq = n(), .groups = "drop")

constraintCodes2 <- constraintCodes[, c("constraintabb", "constraintcategory")]
constraintCodes2 <- head(constraintCodes2, 39) 
constraintCodes2 <- constraintCodes2 %>% 
  pivot_wider(names_from = constraintabb, values_from = constraintcategory) %>%
  pivot_longer(cols = everything(), values_to = "constraintcategory", names_to = "constraintabb") %>%
  mutate(constraintabb = ifelse(constraintabb %in% c("blank1", "blank2", "blank3", 
                                                     "blank4", "blank5", "blank6", 
                                                     "blank7", "blank8"), NA, constraintabb))

levels_panel <- right_join(levels_panel, constraintCodes2, by = "constraintabb") %>%
  mutate(freq = ifelse(is.na(freq), 0, freq))


all_panels <- unique(df$panel)
all_constraints <- constraintCodes2 %>% 
  filter(!is.na(constraintabb)) %>% 
  distinct(constraintabb, constraintcategory)

levels_panel <- expand.grid(
  panel = all_panels, 
  constraintabb = all_constraints$constraintabb,
  stringsAsFactors = FALSE
) %>%
  left_join(all_constraints, by = "constraintabb") %>%
  left_join(levels_panel %>% select(panel, constraintabb, freq), 
            by = c("panel", "constraintabb"))

levels_panel$freq[is.na(levels_panel$freq)] <- 0
levels_panel$panel <- factor(levels_panel$panel, 
                             levels = c(top5_countries, "Other"))


desired_category_order <- c("Economic", 
                            "Governments & Organizations", 
                            "Knowledge",
                            "Material Inputs", 
                            "Negative Side Effects", 
                            "Policies & Rules",
                            "Social & Behavioral")

levels_panel <- levels_panel %>%
  mutate(constraintcategory = factor(constraintcategory, 
                                     levels = desired_category_order))




my_colors <- c(
  "#4daf4a",  # Economic
  "#ffff33",  # Governments & Organizations
  "#984ea3",  # Knowledge
  "#377eb8",  # Material Inputs
  "#ff7f00",  # Negative Side Effects
  "#e41a1c",  # Policies & Rules
  "#a65628"   # Social & Behavioral
)

plot_pie_panel <- function(panel_name) {
  data_panel <- levels_panel %>%
    filter(panel == panel_name, !is.na(constraintcategory)) %>%
    group_by(constraintcategory) %>%
    summarise(total_freq = sum(freq), .groups = "drop")
  
  p <- ggplot(data_panel, aes(x = "", y = total_freq, fill = constraintcategory)) +
    geom_bar(stat = "identity", width = 1) +
    coord_polar("y", start = 0) +
    scale_fill_manual(values = my_colors, na.value = "grey90") +
    theme_void() +
    labs(fill = "Constraint Category", title = panel_name) +
    geom_text(aes(label = total_freq), 
              position = position_stack(vjust = 0.5), 
              size = 4) +
    theme(plot.title = element_text(size = 14, face = "bold"))
  
  return(p)
}


p1 <- plot_pie_panel(top5_countries[1])
p2 <- plot_pie_panel(top5_countries[2])
p3 <- plot_pie_panel(top5_countries[3])
p4 <- plot_pie_panel(top5_countries[4])
p5 <- plot_pie_panel(top5_countries[5])
p6 <- plot_pie_panel("Other")


figs6_pies <- (p1 + p2 + p3) / (p4 + p5 + p6) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

figs6_pies


```




##Figure S7: Estimated CSP affected by each constraint
```{r}
constraint_co2e <- df %>%
  distinct(constraintabb, ISOCode) %>%
  left_join(Naturebase, by = "ISOCode") %>%
  group_by(constraintabb) %>%
  summarise(total_NCS_CO2e = sum(total_ncs_co2e_total, na.rm = TRUE), .groups = "drop")

constraintCodes2 <- constraintCodes[, c("constraintabb", "constraintcategory")]
constraintCodes2 <- head(constraintCodes2, 39)

constraintCodes2 <- constraintCodes2 %>%
  pivot_wider(names_from = constraintabb, values_from = constraintcategory) %>%
  pivot_longer(cols = everything(), values_to = "constraintcategory", names_to = "constraintabb")

constraintCodes2 <- constraintCodes2 %>%
  mutate(constraintabb = ifelse(constraintabb %in% c("blank1", "blank2", "blank3", "blank4", "blank5", "blank6", "blank7", "blank8"), NA, constraintabb))

levels1 <- right_join(constraint_co2e, constraintCodes2, by = "constraintabb")

levels1 <- levels1 %>%
  mutate(total_NCS_CO2e = ifelse(is.na(constraintabb), 0, total_NCS_CO2e))

levels <- levels1 %>%
  arrange(constraintcategory, total_NCS_CO2e)

testLabels <- levels$constraintabb
testLabels <- testLabels[-c(47, 48)]  

levels1$constraintabb <- factor(levels1$constraintabb, levels$constraintabb)

figs7<- levels1 %>%
  filter(!is.na(constraintcategory)) %>%
  filter(!is.na(constraintabb)) %>%
  group_by(constraintabb) %>%
  mutate(
    constraintabb = fct_rev(constraintabb),
    total_NCS_CO2e_Mt = total_NCS_CO2e / 1e6  
  ) %>%
  ggplot(aes(x = constraintabb, y = total_NCS_CO2e_Mt, fill = constraintcategory)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = rev(c(
    "#a65628", "#e41a1c", "#ff7f00", "#377eb8", "#984ea3", "#ffff33", "#4daf4a"
  ))) +
  scale_x_discrete(labels = rev(testLabels)) +
  coord_flip() +
  theme_bw() +
  theme(legend.title = element_text(size = 8)) +
  xlab("Constraint") +
  ylab("Total CMP Potentially Affected by Constraint (Mt CO2e)") +  
  labs(fill = "Constraint Category") +
  geom_text(aes(label = round(total_NCS_CO2e_Mt, 1)), hjust = 1.25, colour = "lightgrey", size = 2.5)

print(figs7)



```


##Figure S8: Jaccard and Euclidean within and between geographies
```{r}

calculate_jaccard_similarity <- function(df, region_col) {
  country_table <- df %>%
    distinct(Country, constraintabb, !!sym(region_col)) %>%
    mutate(presence = 1) %>%
    spread(key = constraintabb, value = presence, fill = 0)

  region_data <- country_table %>%
    select(-Country, -!!sym(region_col))

  regions <- country_table[[region_col]]

  distance_matrix <- vegdist(region_data, method = "jaccard")
  similarity_matrix <- 1 - as.matrix(distance_matrix)

  similarity_df <- as.data.frame(similarity_matrix)
  similarity_df$Country <- rownames(similarity_df)
  similarity_df[[region_col]] <- regions

  similarity_long <- melt(similarity_df, id.vars = c("Country", region_col), variable.name = "Country2", value.name = "Similarity")

  similarity_long[[paste0(region_col, "2")]] <- regions[match(similarity_long$Country2, similarity_df$Country)]

  similarity_long$Comparison <- ifelse(similarity_long[[region_col]] == similarity_long[[paste0(region_col, "2")]],
                                       paste0("Within-", region_col), paste0("Between-", region_col))

  return(similarity_long)
}

calculate_euclidean_distances <- function(df, region_col) {
  country_table <- df %>%
    group_by(Country, constraintabb, !!sym(region_col)) %>%
    summarise(count = n()) %>%
    ungroup()

  total_counts <- country_table %>%
    group_by(Country) %>%
    summarise(total_count = sum(count))

  proportion_table <- country_table %>%
    left_join(total_counts, by = "Country") %>%
    mutate(proportion = count / total_count) %>%
    select(-count, -total_count) %>%
    spread(key = constraintabb, value = proportion, fill = 0)

  constraint_data <- proportion_table %>%
    select(-Country, -!!sym(region_col))

  regions <- proportion_table[[region_col]]

  distance_matrix <- dist(constraint_data, method = "euclidean")

  similarity_df <- as.data.frame(as.matrix(distance_matrix))
  similarity_df$Country <- rownames(similarity_df)
  similarity_df[[region_col]] <- regions

  similarity_long <- melt(similarity_df, id.vars = c("Country", region_col), variable.name = "Country2", value.name = "Distance")

  similarity_long[[paste0(region_col, "2")]] <- regions[match(similarity_long$Country2, similarity_df$Country)]

  similarity_long$Comparison <- ifelse(similarity_long[[region_col]] == similarity_long[[paste0(region_col, "2")]],
                                       paste0("Within-", region_col), paste0("Between-", region_col))

  return(similarity_long)
}

similarity_long_sdgrgion <- calculate_jaccard_similarity(df, "SDGrgion")
similarity_long_subregion <- calculate_jaccard_similarity(df, "Subregion")

similarity_long_sdgrgion_euc <- calculate_euclidean_distances(df, "SDGrgion")
similarity_long_subregion_euc <- calculate_euclidean_distances(df, "Subregion")

p1 <- ggplot(similarity_long_sdgrgion, aes(x = Comparison, y = Similarity)) +
  geom_boxplot() +
  labs(x = "SDG Region Comparison", y = "Jaccard Similarity") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

p2 <- ggplot(similarity_long_subregion, aes(x = Comparison, y = Similarity)) +
  geom_boxplot() +
  labs(x = "UN Subregion Comparison", y = "Jaccard Similarity") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

p3 <- ggplot(similarity_long_sdgrgion_euc, aes(x = Comparison, y = Distance)) +
  geom_boxplot() +
  labs(x = "SDG Region Comparison", y = "Euclidean Distance") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

p4 <- ggplot(similarity_long_subregion_euc, aes(x = Comparison, y = Distance)) +
  geom_boxplot() +
  labs(x = "UN Subregion Comparison", y = "Euclidean Distance") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

grid.arrange(p1, p2, p3, p4, ncol = 2)


```


##Figure S9: Pathway Jaccard similarity

```{r}


df_disaggregated <- rawDF %>% filter(secondscreen == 1, reviewscreen == 0, yearpublication >= 2020) %>% dplyr::select(constraintID, countryconstraint, constraint, constraintpathway, solutions)

df_disaggregated <- df_disaggregated %>% mutate(countryconstraint = strsplit(as.character(countryconstraint), ", ")) %>%  unnest(countryconstraint)

df_disaggregated <- df_disaggregated %>% mutate(constraintpathway = strsplit(as.character(constraintpathway), ", ")) %>%
  unnest(constraintpathway)

colnames(df_disaggregated) <- c("constraintID", "ISOCode", "constraint", "pathway", "solutions")

df_disaggregated <- merge(df_disaggregated, placeCodes , by = "ISOCode")

constraintCodes <- rename(constraintCodes, constraint =constraintcategory)

constraintCodes <- head(constraintCodes, 39)


df_disaggregated <- merge(df_disaggregated,constraintCodes, by = "constraint", all = TRUE)

df_disaggregated <- df_disaggregated %>% filter(pathway != "FPM", pathway != "RWF")


# Create the subregion_pathway column
df_disaggregated <- df_disaggregated %>%
  mutate(subregion_pathway = paste(Subregion, pathway, sep = "_"))

# Create a binary presence/absence table for constraints 
subregion_pathway_table <- df_disaggregated %>%
  distinct(subregion_pathway, constraintabb, Subregion, pathway) %>%
  mutate(presence = 1) %>%
  spread(key = constraintabb, value = presence, fill = 0)

# Function to calculate similarity within each subregion
calculate_similarity_within_subregion <- function(subregion) {
  subregion_data <- subregion_pathway_table %>%
    filter(Subregion == subregion) %>%
    select(-subregion_pathway, -Subregion, -pathway)
  
  subregion_pathway <- subregion_pathway_table %>%
    filter(Subregion == subregion) %>%
    pull(subregion_pathway)
  
  # Calculate Jaccard distance matrix
  distance_matrix <- vegdist(subregion_data, method = "jaccard")
  
  # Convert distance matrix to similarity matrix
  similarity_matrix <- 1 - as.matrix(distance_matrix)
  rownames(similarity_matrix) <- subregion_pathway
  colnames(similarity_matrix) <- subregion_pathway
  
  return(similarity_matrix)
}

subregions <- unique(df_disaggregated$Subregion)

similarity_matrices <- lapply(subregions, calculate_similarity_within_subregion)
names(similarity_matrices) <- subregions

pathways <- unique(df_disaggregated$pathway)
n_pathways <- length(pathways)
sum_similarity_matrix <- matrix(0, nrow = n_pathways, ncol = n_pathways)
count_matrix <- matrix(0, nrow = n_pathways, ncol = n_pathways)
rownames(sum_similarity_matrix) <- pathways
colnames(sum_similarity_matrix) <- pathways
rownames(count_matrix) <- pathways
colnames(count_matrix) <- pathways

# Sum all similarity matrices and count occurrences
for (sim_matrix in similarity_matrices) {
  for (i in 1:nrow(sim_matrix)) {
    for (j in 1:ncol(sim_matrix)) {
      pathway_i <- strsplit(rownames(sim_matrix)[i], "_")[[1]][2]
      pathway_j <- strsplit(colnames(sim_matrix)[j], "_")[[1]][2]
      
      # Debugging: print pathway_i and pathway_j
      print(paste("pathway_i:", pathway_i, "pathway_j:", pathway_j))
      
      if (!is.na(sim_matrix[i, j])) {
        if (pathway_i %in% pathways && pathway_j %in% pathways) {
          sum_similarity_matrix[pathway_i, pathway_j] <- sum_similarity_matrix[pathway_i, pathway_j] + sim_matrix[i, j]
          count_matrix[pathway_i, pathway_j] <- count_matrix[pathway_i, pathway_j] + 1
        } else {
          print(paste("Error: pathway_i or pathway_j not in pathways:", pathway_i, pathway_j))
        }
      }
    }
  }
}

# Debugging: check the final matrices
print("Sum similarity matrix:")
print(sum_similarity_matrix)
print("Count matrix:")
print(count_matrix)


average_similarity_matrix <- sum_similarity_matrix / count_matrix

average_similarity_matrix[is.na(average_similarity_matrix)] <- 0


# Plot heatmap of the average similarity matrix
pheatmap(average_similarity_matrix, 
         cluster_rows = TRUE, 
         cluster_cols = TRUE, 
         display_numbers = TRUE, 
         color = colorRampPalette(c("white", "blue"))(100))

```

##Figure S10: Pathway Euclidean similarity

```{r}

df_disaggregated <- df_disaggregated %>%
  mutate(subregion_pathway = paste(Subregion, pathway, sep = "_"))

df_disaggregated <- df_disaggregated %>%
  group_by(subregion_pathway) %>%
  mutate(total_constraints = n()) %>%
  ungroup() %>%
  group_by(subregion_pathway, constraintabb) %>%
  summarise(proportion = n() / mean(total_constraints)) %>%
  ungroup()

subregion_pathway_table <- df_disaggregated %>%
  spread(key = constraintabb, value = proportion, fill = 0)

subregion_pathway <- subregion_pathway_table$subregion_pathway

# Function to calculate similarity within each subregion using Euclidean distance
calculate_similarity_within_subregion <- function(subregion) {
  subregion_data <- subregion_pathway_table %>%
    filter(grepl(subregion, subregion_pathway)) %>%
    select(-subregion_pathway)
  
  subregion_pathway <- subregion_pathway_table %>%
    filter(grepl(subregion, subregion_pathway)) %>%
    pull(subregion_pathway)
  
  # Calculate Euclidean distance matrix
  distance_matrix <- vegdist(subregion_data, method = "euclidean")
  
  # Convert distance matrix to similarity matrix
  similarity_matrix <- 1 / (1 + as.matrix(distance_matrix))
  rownames(similarity_matrix) <- subregion_pathway
  colnames(similarity_matrix) <- subregion_pathway
  
  return(similarity_matrix)
}

subregions <- unique(gsub("_.*", "", subregion_pathway))

similarity_matrices <- lapply(subregions, calculate_similarity_within_subregion)
names(similarity_matrices) <- subregions

pathways <- c("AWC", "CWR", "AGC", "GrR", "APC", "PeR", "CSF", "AgFo", "RFo", "AFC")
n_pathways <- length(pathways)
sum_similarity_matrix <- matrix(0, nrow = n_pathways, ncol = n_pathways)
count_matrix <- matrix(0, nrow = n_pathways, ncol = n_pathways)
rownames(sum_similarity_matrix) <- pathways
colnames(sum_similarity_matrix) <- pathways
rownames(count_matrix) <- pathways
colnames(count_matrix) <- pathways

for (sim_matrix in similarity_matrices) {
  for (i in 1:nrow(sim_matrix)) {
    for (j in 1:ncol(sim_matrix)) {
      pathway_i <- strsplit(rownames(sim_matrix)[i], "_")[[1]][2]
      pathway_j <- strsplit(colnames(sim_matrix)[j], "_")[[1]][2]
      if (!is.na(sim_matrix[i, j])) {
        sum_similarity_matrix[pathway_i, pathway_j] <- sum_similarity_matrix[pathway_i, pathway_j] + sim_matrix[i, j]
        count_matrix[pathway_i, pathway_j] <- count_matrix[pathway_i, pathway_j] + 1
      }
    }
  }
}

average_similarity_matrix <- sum_similarity_matrix / count_matrix
average_similarity_matrix[is.na(average_similarity_matrix)] <- 0
diag(average_similarity_matrix) <- 1


pheatmap(average_similarity_matrix, 
         cluster_rows = TRUE, 
         cluster_cols = TRUE, 
         display_numbers = TRUE, 
         color = colorRampPalette(c("white", "blue"))(100),
         labels_row = pathways,
         labels_col = pathways)

```





##Figure S11: Sankey diagram 
```{r}

category_colors <- c(
  "Policies & Rules"           = "#e41a1c",
  "Negative Side Effects"      = "#ff7f00",
  "Material Inputs"            = "#377eb8",
  "Knowledge"                  = "#984ea3",
  "Governments & Organizations"= "#ffff33",
  "Economic"                   = "#4daf4a",
  "Social & Behavioral"        = "#a65628")


df_grouped <- df %>%
  group_by(pathway, SDGrgion, constraintcategory) %>%
  summarise(Freq = n(), .groups = "drop")  %>%
  mutate(pathway = as.factor(pathway),
    SDGrgion = as.factor(SDGrgion),
    constraintcategory = as.factor(constraintcategory),
    Freq = as.numeric(Freq))


ggplot(df_grouped, aes(axis1 = pathway, axis2 = SDGrgion, y = Freq)) +
  geom_alluvium(aes(fill = constraintcategory), width = 0.3, alpha = 0.8) +  
  geom_stratum(width = 0.3) +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 3) +
  scale_fill_manual(values = category_colors, na.value = "grey") +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
        axis.title.y = element_text(face = "bold")  
  ) +
  labs(
    x = NULL,
    y = "Frequency",
    fill = "Constraint Category"
  ) +
annotate("text", x = 1, y = -max(df_grouped$Freq) , label = "Pathway", size = 4, fontface = "bold") +
  annotate("text", x = 2, y = -max(df_grouped$Freq) , label = "SDG Region", size = 4, fontface = "bold")




```
